
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ESERCIZIO 2  NOSTRA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%aggiornamento dropbox

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% HIDDEN MARKOV MODEL
            
            % Filtro di viterbi
            % Stima di viterbi
            % Algoritmo di alpha-beta
            % Stringa ricevuta
            % Valutazione delle prestazioni
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%Inizializzazione:

clear all
clc

% Numero delle coppie di bit trasmesse
T=20;

% Combinazioni degli stati nascosti
v1=[0 0];
v2=[0 1];
v3=[1 0];
v4=[1 1];

v=[v1' v2' v3' v4'];

% Osservazioni possibili
w=v;

% Matrice transizione di stato
A=[0.25 0.25 0.5 0; 0.5 0.5 0 0; 0.5 0 0.25 0.25; 0.25 0 0 0.75];

% Probabilità dello stato iniziale
P0=[0.25 0.25 0.25 0.25]';

% Rumore del canale
epsilon = [0.05 0.25 0.45];

% Ciclo al variare del rumore epsilon di canale
for u=1:3

err_seq_fil_viterbi1 =0;
err_stima_viterbi1 =0;
err_alg_alfabeta1=0; 
err_stringa1=0;

num_prob_non_0=zeros(1,4);
coeff_z=zeros(1,4);
    
% Matrice transizione stato/uscita
% Calcolo della distanza di Hamming e della matrice C
for l=1:4
   for k=1:4
      d=haming(v(:,l),w(:,k),1);
      C(l,k)=(1-epsilon(u))^(2-d)*epsilon(u)^d; 
   end
end

% Generazione di una sequenza di Markov e delle sue osservazioni
[oss,stat] = hmmgenerate(T,A,C); %funzione matlab

% Ciclo per effettuare 100 realizzazioni dei 4 stimatori 
for n = 1:100

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo del filtro alla viterbi: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

seq_fil_vit = [];

for t=1:T
    
seq_fil_vit = [seq_fil_vit, viterbone(t,P0,oss,A,C)'];
[maxProb(t) q_ott(t)]=max(viterbone(t,P0,oss,A,C));

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stima di viterbi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Si calcola la stima di viterbi a partire dai dati ottenuti dal filtro
% di viterbi. L'algoritmo parte da T e si muove all'indietro


prob_stimaVit = maxProb;
oss_stima = q_ott;
for t=T:-1:2
 
    for i=1:4
        vett_stima_temp(i,t) = seq_fil_vit(i,t-1)*A(i,oss_stima(t));
    end
    [prob_stimaVit(t-1) oss_stima(t-1)] = max(vett_stima_temp(:,t));
    
end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Algoritmo alpha-beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global betaTemp;
betaTemp=zeros(4,T);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Passata all'indietro: Beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for t=T:-1:1
if t==T
     for i=1:4
        betaTemp(i,t)=1;
     end
  else
     for i=1:4
         for j=1:4
         betaTemp(i,t)=betaTemp(i,t)+ (A(i,j)*C(j,oss(t+1))*betaTemp(j,t+1));
         end
      end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Passata in avanti: Alpha
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alpha
alpha=zeros(4,T);

for i=1:4
           alpha(i,1)=P0(i)*C(i,oss(1)); 
end 

for t=1:T-1
    for j=1:4 
        for i=1:4
          alpha(j,t+1) = alpha(j,t+1)+(A(i,j)*C(j,oss(t+1))*alpha(i,t));
        end
    end    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fondo le due passate e la chiamo alfa_beta
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alfa_beta=zeros(4,T);

for t=1:T
   for i=1:4
      den=0;
      num=alpha(i,t)*betaTemp(i,t);
      for j=1:4
         den=den + alpha(j,t)*betaTemp(j,t);
      end
      alfa_beta(i,t)=num/den;
   end
end

% ricerca della sequenza ottima
for k=1:T

[alg_beta(k) alg_alfa(k)]=max(alfa_beta(:,k));

end

%%%%%%%%%%%%%%%%%%%%%%%
%Stringa di osservazioni
%%%%%%%%%%%%%%%%%%%%%%%

stringa=oss;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verifica delle prestazioni:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Riconversione dei stati in bit:

stat1 = conv_bit(stat,T);
q_ott1 = conv_bit(q_ott,T);
oss_stima1 = conv_bit(oss_stima,T);
alg_alfa1 =  conv_bit(alg_alfa,T);
oss1 = conv_bit(oss,T);



% Bit error rate dopo la stima estimation

% Calcolo la distanza di Hamming tra gli stati e gli stimatori

d_ott = haming(stat1,q_ott1,T);
d_stima = haming(stat1,oss_stima1,T);
d_alfa = haming(stat1,alg_alfa1,T);
d_oss = haming(stat1,oss1,T);

d_ott1 = 0;
d_stima1 = 0;
d_alfa1 = 0;
d_oss1 = 0;

% Sommo le distanze di Hamming per tutti i T-istanti

for i =1:T
    d_ott1 = d_ott1 + d_ott(i);
    d_stima1 = d_stima1 + d_stima(i);
    d_alfa1 = d_alfa1 + d_alfa(i);
    d_oss1 = d_oss1 + d_oss(i);    
end

% Calcolo bit error rate dividendo la distanza di Hamming per 2*T

err_seq_fil_viterbi(n)   = d_ott1/(2*T);
err_stima_viterbi(n)   = d_stima1/(2*T);
err_alg_alfabeta(n)  = d_alfa1/(2*T);
err_stringa(n) =  d_oss1/(2*T);

% Sommo la bit error rate per tutte le 100 realizzazioni

err_seq_fil_viterbi1   = err_seq_fil_viterbi1 + err_seq_fil_viterbi(n);
err_stima_viterbi1   = err_stima_viterbi1 + err_stima_viterbi(n);
err_alg_alfabeta1  = err_alg_alfabeta1 + err_alg_alfabeta(n);
err_stringa1 =  err_stringa1 + err_stringa(n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calcolo della probabilit� della sequenza stimata Log della P(Xi=Yi|Zi)_
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 log_p=P0(1)*ones(1,4);

for t=2:T-1
    %Probabilit�delle sequenze dei 4 stimatori
    
    % Prob. di partire dall' istante iniziale e passare attraverso 
    % la sequenza pi� probabile verso l'istante finale
    
   log_p(1)=log_p(1)*A(q_ott(t),q_ott(t+1))*C(q_ott(t),oss(t));
   log_p(2)=log_p(2)*A(oss_stima(t),oss_stima(t+1))*C(oss_stima(t),oss(t));
   log_p(3)=log_p(3)*A(alg_alfa(t),alg_alfa(t+1))*C(alg_alfa(t),oss(t));
   log_p(4)=log_p(4)*A(oss(t),oss(t+1))*C(oss(t),oss(t));
end

% Calcolo del logaritmo
for i=1:4
  log_p(i)=log(log_p(i));
end

% Conto il numero di elementi di probabilit� piccolissima

 for i=1:4
     if log_p(i) ~= -Inf
       num_prob_non_0(i)=num_prob_non_0(i)+1;
       coeff_z(i)=coeff_z(i) + log_p(i);
     end
   end

  
end


for i=1:4
  coeff_z(i)=coeff_z(i)/num_prob_non_0(i);
end


disp('Indici con epsilon:')

disp(epsilon(u)) 

disp('Bit error rate:')
disp('filtro - stimatore - alpha-beta - Stringa')
disp([err_seq_fil_viterbi1/100 err_stima_viterbi1/100 err_alg_alfabeta1/100 err_stringa1/100])

disp('Probabilit� della sequenza stimata:')
disp('filtro - stimatore - alpha-beta - Stringa')
disp([coeff_z(1) coeff_z(2) coeff_z(3) coeff_z(4)])

disp('Stati Filtro Stimatore Alfa-Beta Osservazioni')
disp([stat' q_ott' oss_stima' alg_alfa' oss'])

end
